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Abstract 

We present a method to measure cosmic magnetic fields with ultra high energy 
cosmic rays (UHECRs) . We apply an advanced autocorrelation method to sim- 
ulated UHECRs which includes their directional as well as energy information. 
Without explicit knowledge of the UHECR sources, such measurements are sen- 
sitive to the number of sources and to the magnetic field strength subjected to 
the UHECRs. Using a UHECR Monte Carlo model including sources, random 
walk propagation and a coherent deflection, we explain the procedure of recon- 
structing the allowed phase space of the model parameters from a simulated 
autocorrelation measurement. 
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1. Introduction 

Recent results of the Pierre Auger Observatory P, Q show a correlation 
of ultra high energy cosmic rays (UHECRs) above 56 EeV with nearby active 
galactic nuclei (AGN) at distances up to 75 Mpc. The exact origin of the 
UHECRs and the positions of their sources still remain unknown but it is favored 
that UHECRs are accelerated at discrete sources. 

The UHECRs propagate through extragalactic and galactic magnetic fields, 
and consequently carry information about the direction and strength of the 
fields. For UHECRs originating from a single source with energies above a 
certain threshold, the influence of the fields is expected to result in an energy 
ordering with respect to the direction of the source. 

Galactic magnetic fields are known to some extent, from different measure- 
ment methods. For recent reviews refer to, e.g. js, 4]. Information about extra- 
galactic fields exists for a number of galaxy clusters, however, the knowledge of 
magnetic fields in filaments is uncertain to at least three orders of magnitude 
[1, @, 01 ■ Measurements of UHECRs therefore have good potential of providing 
additional information on these magnetic fields. 
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In this contribution we present a new method of obtaining information on 
the magnetic fields without prior knowledge of the source positions. This is 
achieved by a statistical approach based on the concept of energy dependent 
angular ordering of the UHECRs mentioned above. This publication is orga- 
nized as follows. First we define an energy-energy-correlation observable used 
to measure the angular ordering. Then we explain the analysis procedure, and 
demonstrate its application using a Monte Carlo generated UHECR scenario. 
Finally, we present a method to test models of UHECR emission and propaga- 
tion using a simulated energy-energy-correlation measurement. In this context 
we reconstruct the parameters of a simple magnetic field model which are the 
number of sources emitting UHECRs, and the turbulent magnetic field strength, 
while keeping a coherent magnetic field strength constant. 



2. Definition of Energy-Energy-Correlations 

Energy-energy-correlations are a well known quantity from high energy phys- 
ics, see e.g. Their definition usually includes the energies and the angular 
distances of particles, which enables investigation of energy ordering in the 
UHECR sky. We define the energy-energy-correlation ftij between the UHECRs 
i and j by 



ma,) - { E{a,))){E,{a, ) - {E{a,))) 

Here Ei is the energy of the UHECR z, and ai denotes its angular distance 
with respect to the center of a region of interest (ROI) (figure [Ij . The ROI 
covers a limited solid angle and will be precisely outlined in section [3^ {E{ai)) 
denotes the corresponding average energy value of the UHECRs arriving within 
the same ring interval relative to the ROI center. 

Using ilij, angular ordering of the UHECRs is measured in the following 
sense, 

1. A pair of UHECRs, one being above and the other below the corresponding 
average energy values, results in a negative correlation flij < 0. This is a 
typical case for a background contribution. 

2. A pair with both UHECRs having energy values above or below the aver- 
age energy value at their corresponding angular distance gives a positive 
correlation ilij > 0. Here signal UHECRs are expected to contribute. 

This means we expect overall larger flij for a mixing of coherently arriving 
UHECRs from a source with background UHECRs, than for an exclusively 
isotropic arrival. In this way, the energy-energy-correlations as defined above 
quantify the coherence in angular ordering of the UHECRs coming from the 
same source, and provide a separation of incoherently arriving UHECRs from 
different sources on a statistical basis. 
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Figure 1: Schematic view of energy-energy-correlations 

3. Analysis Method 

3.1. Signal Data Set 

In order to evaluate the sensitivity of the energy-energy-correlation observ- 
able, we use instead of real data, a set of simulated UHECRs resulting from a 
source distribution. We will reference this in the following as signal data set. 

Future Monte Carlo simulations of individual UHECRs aim to include spe- 
cific candidate sources, a realistic distribution of galaxies, galaxy clusters and 
filaments, and a proper propagation through the corresponding media and mag- 
netic fields. Here we use a simple representation of such a UHECR Monte Carlo 
generator which includes a random walk propagation and a coherent deflection 
of the UHECRs. Random walk propagation can be considered to approximate 
the influence of extragalactic mag netic fields (EGMF) [§], while it is assumed 
that the dominating effect of the galactic magnetic field (GMF) is a coherent 
deflection [lo| . 

The random walk propagation is simulated by varying the angle 9 between 
the source direction and the arrival direction of the UHECR. depends on the 
distance D which a proton propagates through magnetic fields of strength B in 
zones of the correlation length A. The angle is chosen randomly according to a 
2-dimensional Gaussian distribution with the width 

- 0.025 deg (f ) ^'^ (^) (^) (^) ' . (2) 

For the GMF we adopt a homogeneous magnetic field perpendicular to the 
galactic plane as a simple representative model which simulates its most impor- 
tant features. In galactic coordinates this results in the following deflection. 
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Ai" ~ CcoherentField 

Every UHECR is shifted in galactic longitude i depending on its energy by 
the amount A^. The galactic latitude remains unchanged. Following from this 
the coherent deflection strength is maximal near the galactic plane, while it 
drops to zero at the poles. Although this gives a completely difii'erent global 
deflection pattern of UHECRs compared to e.g. [13], the local effect of an 
energy ordering is the same. For the signal data set as well as for all other 
Monte Carlo sets used in this paper we adopt a coherent deflection strength 
of CcoherentField — 10 rad. In casc of an UHECR of 60EeV this corresponds 
to a maximal deflection of 10°, and an average deflection of 7.5° which is 
compatible with the results of jlO| . 

For the energy distribution of the UHECRs we follow the spectrum measured 
by the Pierre Auger Observatory [lH, 



J{E) oc (4) 

with 7 = -2.69 for 4 EeV < S < 40 EeV, 7 = -4.2 for S > 40 EeV. 

In the production of our signal data set wc first generated 10 randomly 
distributed sources in the sky. We assumed that all sources have a similar 
distance to the earth such that equation [5] simplifies to 

ffl — CfiandomField ( i nlS _ V 

(5) 

We produced the signal data set with C RandomFieid = 10 rad. A scenario 
compatible with this choice is, e.g. D = 50 Mpc, A = 1 Mpc, and B = ?> nC 
We further assumed that all sources have identical luminosity, and each of them 
emitted 1000 protons. The total signal data set amounts to 10000 UHECRs 
with energies above 5 EeV. This energy threshold has been chosen to be above 
the ankle measured in [lli, such that the UHECRs can be considered to have 



extragalactic origin [Ij]. In figure 2(a) we show the positions of the sources 



(asterisk symbols), and the resulting UHECR distribution (point symbols) as 
observed at the earth using the Hammer projection of the galactic coordinates. 
The UHECR distribution appears to be almost isotropic. 



3.2. Regions oj Interest 

We define a ROI as a region that is close to a source candidate of UHECRs. 
Since no UHECR source has been clearly identified so far, we applied a simple 
iterative cone algorithm to find ROFs that have an increased probability to 
contain a source. The algorithm works as follows. 

1. Select all UHECRs with energies above E^i^ = 60 EeV as initial seeds, 

2. For all seeds, define a corresponding ROI by assigning all UHECRs with 
angular distances less than amax = 0-2 rad. 
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Figure 2: |(a)| Arrival directions of the signal data set (black point symbols) and the sources 
(red asterisk symbols) , |(b)| UHECRs belo nging to regions of interest (black point symbols) 
and the sources (red asterisk symbols) , | (c) [ Reconstruction quality of the source direction using 
the cone algorithm (red line: 68%-quantile) , f(d) | Energy-energy-correlations of the signal data 
set 
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3. Calculate the center of mass of each ROI using the energies of the UHECRs 
as weights, 

4. Use the center of mass of the ROI as a new seed, and iterate starting from 
item 2. 

Note that in this algorithm every UHECR can be part of several ROIs. The 
algorithm is processed in total three times, and the last resulting ROIs are taken 
for further analysis. 

In order to test the accuracy of reconstruction of a source direction with the 
ROI method, we produced 100 additional simulated data sets using the same 
magnetic field strength, and the same number of sources as above. We varied 
the random generator seeds to produce different source positions and UHECRs. 



The performance of the ROI algorithm is shown in figures 2(b) and 2(c) In 
figure [2(b)] one can see, that in the signal data set 9 out of 10 sources are covered 
by ROIs. Using the 100 additional data sets with the identical parameters, 95% 



of the sources are covered by ROIs. In figure 2(c) the angular distance Aa 
between a source and the closest ROI is shown. We obtain an angular resolution 
in terms of the 68%-quantile of 0.12 rad. Note that some reconstructed ROIs 
do not correspond to sources, and provide a background contribution to the 
following analysis. 

3.3. Calculation of Energy- Energy- Correlations 

In the next step we calculate the energy-energy-correlations of all pairs of 
UHECRs that belong to the same ROI. Each value of flij (equation [1]) is filled 
into a histogram at both angular distances and aj . Finally, in every angular 
bin of a we calculate the average value and its uncertainty. In figure 2(d) we 



show the distribution of fl where we have included the flij values of all ROIs 
simultaneously. The corresponding error bars are small compared to the size 
of the symbols. The result is the angular distribution of the energy-energy- 
correlations, which represents an estimator of coherence of the UHECR sky. 

In experiments measuring UHECRs, this distribution can be obtained di- 
rectly from the data without any assumption on sources, propagation, or mag- 
netic fields. 



4. Evaluation of UHECR models 

Each state-of-the-art model of UHECR emission and propagation can be 
confronted with the distribution of an energy-energy-correlation measurement 
as presented above. In this section we present a method of evaluating such a 
model using the data distribution. To explain the evaluation procedure we use 
simplified models of the UHECR sky. 

As an example, we demonstrate the potential to exclude a model of isotropic 
arrival directions by comparison with the signal data distribution shown in fig- 
ure 2(d) Furthermore, we evaluate variants of the above mentioned magnetic 
field model using different numbers of sources and magnetic field parameters. 
We demonstrate how phase space regions of possible parameter settings of this 
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Figure 3: |(a)| Energy-energy-correlations of the signal data (black square symbols) compared 
with the r esult of isotropic arriv al di rections (red triangle symbols), [(b)] ^-distribution of the 
first bin of |(a)| (same color code) , [(c)] Error contours for the reconstructed parameters resulting 
from the negative log-likelihood analysis (dark blue = Itr-contour, light blue = 2(T-contour, 
dark green = Scr-contour, light green = 4(7-contour, yellow = 5iT-contour, red > 5<t), the black 
point symbol shows the parameters of the signal data set, |(d)| Error contours using 100 signal 
data sets (same color code as in|(c)| 



model can be excluded. Finally we quantify the precision of the reconstruction 
of these parameters. 

4.1. Isotropic UHECR Arrival Directions 

We compare the energy-energy-correlations of the signal data with the hy- 
pothesis of isotropic arrival. We simulated 100 UHECR data sets with isotropic 
arrival directions. For each set we calculated the energy-energy-correlations 
as described in section |3l The average value of the resulting fl-distributions 



is shown in figure 3(a) as the triangle symbols. The error bars represent the 
spread of the average values of the individual distributions. The square symbols 
show the same signal data distribution as in figure |2(d)[ 

In order to quantify the agreement of the isotropic arrival model with the 
signal data distribution we use the negative log-likelihood method. For the 
comparison of the distributions we take bin- by-bin correlations into account. 
We make use of the approximately Gaussian shape of the VL distributions within 



7 



the angular a bins (e.g. figure [3(b)] ) , and define for tlie negative log-likeliliood 
value 

L = -2\n (^ML^exp (^-^ix {y)f ■ F'^ • {x (y)))) . (6) 
Here x and (y) are n-dimensional vectors (n=10, the number of bins in 



figure 3(a) ) containing the signal data values of each a bin (x), and the mean 



of the isotropic distributions ((y)). V is the covariance matrix of the histogram 



resulting from the isotropic model (figure 3(a)). The coefficients of V of the 
distribution y are defined as follows, 

Vij ^ {ViVj) ~ {y,){yj). (7) 

The resulting value oi L — 295 needs to be compared with an isotropic 
reference distribution in L. Using isotropic distributions instead of the signal 
data we obtain a narrow distribution centered at Lo = —46 with a root mean 
square of i_Rj\/s = 5. These Numbers suggest that an isotropic UHECR scenario 
is excluded. In the analysis below we explicitly show that isotropy is excluded 



by more than five standard deviations (figure 3(c) I 



4.2. Magnetic Field Model 

In this section we confront a large number of variants of the magnetic model 
with the signal data distribution (figure |2(d)| . We vary the turbulent mag- 
netic field as represented by the parameter CRandomFieid of equation [5l and the 
number of sources Nsource emitting UHECRs with the same luminosity. The 
parameter for the regular magnetic field is kept fixed, as we found that the 
presented method has only minor sensitivity to variations of CcoherentFieid at a 
fixed value of Nsource- Wc kept the total number of UHECRs Njjhecr = 10000 
constant such that each source emits NjjHECR/Nsource UHECRs. 

In order to cover the relevant regions of the parameter space we chose an 
almost logarithmic binning with 25 different numbers of sources (between 1 and 
10000), and 16 different values of CaandomFieid (between 1 rad and 100 rad). 
For each of these values we simulated 100 data sets including the random walk 
and the coherent defiection, and determined the corresponding energy-energy- 
correlation distributions. The mean values (y) and the covariance matrix V are 
obtained for each pair of the model parameter values {Nsource, CfiandomFieid)- 

For these values we calculated the negative log-likelihood values 
L{Nsource-, C RandomFieid) accordiug to equation [S] using for x again the sig- 



nal data of figure 2(d) In order to reduce fluctuations resulting from the finite 
number of data sets in each {Nsource, C RandomFieid) bin, we smoothed the bins 
of the L{N source, C RandomFieid) distribution with a 5 X 5-kernel. The center 
weight of the kernel is K^-^ = 5. All neighbor weights of K^-^ have the weight 
K22 = K23 = if 24 = K32 = if 34 = if 42 = if 43 = if 44 = 2, and all second 
neighbors have zero weight except if 13 — if 31 — if 35 = if 53 = 1. This means 
the value of every bin is replaced by a weighted mean of the neighboring bins, 
according to the weights of the kernel. 
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The result of this approximation is shown in figure 3(c) The red area in- 



dicates the region of the model parameters which are excluded at the level of 
more than five standard deviations. The evaluation of this model using the sig- 
nal data distribution together with the log-likelihood method therefore results 
in clear constraints on the allowed phase space of the model parameters. The 
situation of an isotropic sky can be found at N source — 10000 sources which is 
well within the excluded phase space region. 

Beyond this, we have studied the reconstruction quality of the most likely pa- 
rameter settings of the model under evaluation. As the reference, the true input 



parameters of the signal data set are shown by the point symbol in figure 3(c) 



The contours in figure 3(c) represent levels of n standard deviations where 
the inner dark blue region gives the result for n = 1 standard deviation. The 
input parameter values are reconstructed within a region of two standard de- 
viations. The shape of the distribution exhibits a continuous minimum region 
within the three standard deviation level. 

Above the three standard deviation level, the contour features the following 
reduced sensitivity. For the case of a small number of sources Ngource < 2 an 
enhanced strength of the random field component is able to partly simulate 
effects originating from the Ngource = 10 source scenario. The region of small 
turbulent fields cannot be excluded for the region around the original num- 
ber of sources Ngource ~ 10. Here the coherent component of the field with 
CcoherentFieid = 10 rad gives a dominating effect. 

For evaluating a given model using a measured data distribution of energy- 
energy-correlations, the distribution shown in figure 3(c) and its contours will 
be the main result. 

For further crosschecks of the presented method, we produced 100 additional 
data sets with the same input parameters as the signal data set but with different 
random seeds. For each of these data sets the reconstruction of the parameters 
is repeated as described above. 

In order to test the validity of the error contours we compared the recon- 
structed parameters with the input parameters of the model. We reconstructed 
the input parameter values for 40 data sets within one standard deviation, for 
81 data sets within 2(7, and for all 100 sets within 3(7. These values are consis- 
tent with the expectation for a simultaneous determination of two parameters 
following a Gaussian distribution (39.3% within Icr, 86.5% within 2(7, and 98.9% 
within 3(7). 

We then averaged over the 100 uncertainty contours where the resulting dis- 
tribution is shown in figure |3(d)[ To investigate the potential of a systematic 
shift of the reconstructed parameters introduced by our method we checked the 
minimum of figure |3(d)| which lies within 1 standard deviation of the input pa- 
rameters. This demonstrates that a possible systematic shift is small compared 
to the statistical precision of our procedure. 

In addition we controlled the potential bias introduced by the specific choice 
of the random numbers we used to produce the initial signal data set. A simple 
example for such a bias would be two sources in the same line of sight, which 
would appear as one. Both figures 3(c) and |3(d)| exhibit similar shapes of the 
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allowed parameter space and demonstrate that our choice of the initial data set 
represents a typical realization of a UHECR sky with Ngource = 10 sources, and 
field strength of CRandomFieid = CcoherentFieid = 10 rad. 

5. Conclusions 

In this paper we introduced an energy-energy-correlation observable VL de- 
signed to obtain information on sources and propagation of UHECRs. The Vt- 
distribution can be obtained directly from measured data without assumptions 
on the underlying UHECR sky. 

As models of UHECR emission and propagation become available they can 
be confronted with the O measurements. For quantitative evaluation of the 
models, we have presented a likelihood method. 

To check the sensitivity of our procedure, we compared simple models with a 
simulated 51 measurement. We demonstrated how the phase space of model pa- 
rameters, here number of sources and a parameter reflecting the average random 
magnetic field strength, were constrained using a large number of simulations 
where these parameters were varied. 
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